Plantear el problema de minimos cuadrados como un problema de minimizacion:
Tenemos que \(Y=X\beta + \epsilon\), por lo que donde \(\epsilon\) es el error, por lo que podemos plantear el problema como la minimizaci´ón de este error \[min_{ } ||\epsilon||^2 \] que es equivalente a \[min_{\beta} ||Y-X\beta||^2\] Para encontrar el estimador \(\hat{\beta}\) es necesario derivar e igualar a 0: Derivando respecto a \(\beta\): \[-2X^t (Y-X\hat{\beta})=0\] \[2X^tY=2X^tX\hat{\beta}\] \[\hat{\beta}=(X^tX)^{-1}X^tY\] Es el \(\hat{\beta}\) que soluciona el problema de minimización.
Con esto podemos ver que para cada: \[\hat{\beta}_{i}=\frac{\sum_{i}x_{i}y_{i}}{\sum_{i}x_{i}^2}\] La cuál es una composición lineal de los parametros. Desde el planteamiento estamos buscando que el resultado sea lineal en los parametros, no en los datos por lo que este metodo serviria para poder aproximar polinomios de la orma \(Y=X^2\)
Recordemos que al encontrar la proyección de un vector \(Proy_{X}(y)= \frac{<x,y>}{||x||^2}y\) el error de proyeccion \(Y- Proy_{x}(Y)\) será ortogonal a la proyección. En el caso del problema de minimizar el error, nosotros estamos proyectando la variable dependiente Y sobre el espacio formado por nuestros datos X. y buscamos las combinaciones sobre X que hagan este error mas pequeño.
Viendolo como el teorema de Pitagoras el cu´ál nos dice que el cuadrado de la hipotenisa de un triangulo es igual a la suma de los otros dos lados \(h^2=x^2+y^2\) si el angulo formado entre \(x\) y \(y\) es de 90 grados, en este caso \(H=Y\) las variables dependientes, \(x=Proy_{x}(Y)\) y \(y= error\), entonces la manera en la que se puede cumplir el teorema es que el error sea ortogonal a la proyeccion (angulo de 90 grados) sino el error seráa mas grande.
Al definir una columna de unos estamos cambiando un poco el problema. Antes de agregar la columna: \[y_{j}=\beta_{1}x_{1j}+...+\beta_{n}x_{nj}\] Al definir la columa de unos estamos dando cierta holgura al modelo ya que al agregar un parametro \(\alpha\) no estamos forzando a que la aproximacion pase por el origen. \[y_{j}=\alpha+\beta_{1}x_{1j}+...+\beta_{n}x_{nj}\]
Planteando el problema como \[Y= X\beta +\epsilon\] con \[\epsilon \sim N(0,\sigma^2I_{n})\]
podemos probar que \[Y \sim N(X\beta,\sigma^2I_{n})\]
Demostracion: \[Y= X\beta +\epsilon\] \[\epsilon=Y-X\beta\]
Ahora, falta demostrar que Y es Normal: Por propiedades de la Normal como \(\epsilon \sim N(0,\sigma^2I_{n})\) \[\epsilon + X\beta \sim N(X\beta, \sigma^2I_{n})\] (Esto era m´ás rapido pero bueno…)
Entonces con esto podemos escribir la funcion de verosimilitud de Y como: \[ L(\beta, \sigma^2)=\prod \frac{1}{\sqrt{2\pi\sigma^2}} exp(-\frac{(y-X\beta)^2}{2\sigma^2})\]
\[=(\frac {1}{\sqrt{2\pi\sigma^2}})^n exp (- \frac{1}{2\sigma} (y-X\beta)^t(y-X\beta))\]
A la cual le sacamos Logaritmo para despu´és poder buscar los parametros que maximicen la funcion: \[Log(L(\beta, \sigma^2))=-\frac{n}{2}ln(2\pi)-\frac{n}{2}ln(\sigma^2)-\frac{1}{2\sigma^2}(y-X\beta)^t(y-X\beta)\] Ahora hay que maximizar esa funcion sobre \(\beta\) y \(\sigma^2\) Respecto de \(\beta\) \[\frac{dL}{d\beta}=\frac{1}{\sigma^2}(y-X\beta)^tX=0\] \[Xy=X^tX\beta\] \[\hat{\beta}=(X^tX)^{-1}Xy\]
Respecto de \(\sigma^2\) \[\frac{dL}{d\sigma^2}= - \frac{n}{2\sigma^2}+\frac{1}{2\sigma^4}(y-X\beta)^t(y-X\beta)=0\] \[\sigma^2=\frac{1}{n} (y-X\beta)^t(y-X\beta)\] que puede resolverse usando \(\hat{\beta}\) que encontramos en la solucion anterior.
Como pudimos ver el resultado de maxima verosimilitud \[\hat{\beta}=(X^tX)^{-1}Xy\] Es igual al resultado de \(\hat{\beta}\) de minimos cuadrados
El teorema de Gauss Markov plantea que un modelo de regresion lineal en el que se cumple que:
Entonces el mejor estimador lineal insengado de los oeficientes \(\beta\) es el resultante del problema de minimizar los errores al cuadrado. si existe.
summary(reg_precio)
Call:
lm(formula = precio ~ carat + depth + table + x + y + z, data = diamonds_num)
Residuals:
Min 1Q Median 3Q Max
-23878.2 -615.0 -50.7 347.9 12759.2
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 20849.316 447.562 46.584 < 2e-16 ***
carat 10686.309 63.201 169.085 < 2e-16 ***
depth -203.154 5.504 -36.910 < 2e-16 ***
table -102.446 3.084 -33.216 < 2e-16 ***
x -1315.668 43.070 -30.547 < 2e-16 ***
y 66.322 25.523 2.599 0.00937 **
z 41.628 44.305 0.940 0.34744
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 1497 on 53933 degrees of freedom
Multiple R-squared: 0.8592, Adjusted R-squared: 0.8592
F-statistic: 5.486e+04 on 6 and 53933 DF, p-value: < 2.2e-16
Que tan bien ajusta el modelo depende de la correlacion que tenga cada una de las variables explicativas sobre el precio, esto lo podemos ver a trav´és de:
cor(diamonds_num)
price carat depth table x y z
price 1.0000000 0.92159130 -0.01064740 0.1271339 0.88443516 0.86542090 0.86124944
carat 0.9215913 1.00000000 0.02822431 0.1816175 0.97509423 0.95172220 0.95338738
depth -0.0106474 0.02822431 1.00000000 -0.2957785 -0.02528925 -0.02934067 0.09492388
table 0.1271339 0.18161755 -0.29577852 1.0000000 0.19534428 0.18376015 0.15092869
x 0.8844352 0.97509423 -0.02528925 0.1953443 1.00000000 0.97470148 0.97077180
y 0.8654209 0.95172220 -0.02934067 0.1837601 0.97470148 1.00000000 0.95200572
z 0.8612494 0.95338738 0.09492388 0.1509287 0.97077180 0.95200572 1.00000000
Como podemos ver en la tabla de correlaciones, precio est´á correlacionada altamente con “carat”, “x”,“y”,“z” y tienen cierta correlacion con death y table, esto me dice que las variables escogidas sirven para explicar el precio. Podemos ver esta relacion de manera grafica a continuacion en el siguiente Diagrama de dispersion :
library(GGally)
library(ggplot2)
my_fn <- function(data, mapping, ...){
p <- ggplot(data = data, mapping = mapping) +
geom_point() +
geom_smooth(method=loess, fill="red", color="red", ...) +
geom_smooth(method=lm, fill="blue", color="blue", ...)
p
}
g = ggpairs(diamonds_num,columns = 1:7, lower = list(continuous = my_fn))
g
plot: [1,1] [==-------------------------------------------------------------------------------------------------------] 2% est: 0s
plot: [1,2] [====-----------------------------------------------------------------------------------------------------] 4% est: 4s
plot: [1,3] [======---------------------------------------------------------------------------------------------------] 6% est: 4s
plot: [1,4] [=========------------------------------------------------------------------------------------------------] 8% est: 4s
plot: [1,5] [===========----------------------------------------------------------------------------------------------] 10% est: 4s
plot: [1,6] [=============--------------------------------------------------------------------------------------------] 12% est: 4s
plot: [1,7] [===============------------------------------------------------------------------------------------------] 14% est: 4s
plot: [2,1] [=================----------------------------------------------------------------------------------------] 16% est: 4s
plot: [2,2] [===================--------------------------------------------------------------------------------------] 18% est: 2m
plot: [2,3] [=====================------------------------------------------------------------------------------------] 20% est: 2m
plot: [2,4] [========================---------------------------------------------------------------------------------] 22% est: 2m
plot: [2,5] [==========================-------------------------------------------------------------------------------] 24% est: 1m
plot: [2,6] [============================-----------------------------------------------------------------------------] 27% est: 1m
plot: [2,7] [==============================---------------------------------------------------------------------------] 29% est: 1m
plot: [3,1] [================================-------------------------------------------------------------------------] 31% est: 1m
plot: [3,2] [==================================-----------------------------------------------------------------------] 33% est: 2m
plot: [3,3] [====================================---------------------------------------------------------------------] 35% est: 3m
plot: [3,4] [=======================================------------------------------------------------------------------] 37% est: 2m
plot: [3,5] [=========================================----------------------------------------------------------------] 39% est: 2m
plot: [3,6] [===========================================--------------------------------------------------------------] 41% est: 2m
plot: [3,7] [=============================================------------------------------------------------------------] 43% est: 2m
plot: [4,1] [===============================================----------------------------------------------------------] 45% est: 2m
plot: [4,2] [=================================================--------------------------------------------------------] 47% est: 2m
plot: [4,3] [===================================================------------------------------------------------------] 49% est: 2m
plot: [4,4] [======================================================---------------------------------------------------] 51% est: 3m
plot: [4,5] [========================================================-------------------------------------------------] 53% est: 2m
plot: [4,6] [==========================================================-----------------------------------------------] 55% est: 2m
plot: [4,7] [============================================================---------------------------------------------] 57% est: 2m
plot: [5,1] [==============================================================-------------------------------------------] 59% est: 2m
plot: [5,2] [================================================================-----------------------------------------] 61% est: 2m
plot: [5,3] [==================================================================---------------------------------------] 63% est: 2m
plot: [5,4] [=====================================================================------------------------------------] 65% est: 2m
plot: [5,5] [=======================================================================----------------------------------] 67% est: 2m
plot: [5,6] [=========================================================================--------------------------------] 69% est: 2m
plot: [5,7] [===========================================================================------------------------------] 71% est: 2m
plot: [6,1] [=============================================================================----------------------------] 73% est: 2m
plot: [6,2] [===============================================================================--------------------------] 76% est: 2m
plot: [6,3] [=================================================================================------------------------] 78% est: 2m
plot: [6,4] [====================================================================================---------------------] 80% est: 2m
plot: [6,5] [======================================================================================-------------------] 82% est: 1m
plot: [6,6] [========================================================================================-----------------] 84% est: 1m
plot: [6,7] [==========================================================================================---------------] 86% est: 1m
plot: [7,1] [============================================================================================-------------] 88% est: 1m
plot: [7,2] [==============================================================================================-----------] 90% est:50s
plot: [7,3] [================================================================================================---------] 92% est:41s
plot: [7,4] [===================================================================================================------] 94% est:32s
plot: [7,5] [=====================================================================================================----] 96% est:22s
plot: [7,6] [=======================================================================================================--] 98% est:11s
Graficamente así se ven los valores del modelo contra los observados
data<-as.data.frame(cbind(reg_precio$fitted.values,diamonds_num$price))
p <- plot_ly(data = data, x = ~data[,1], y = ~data[,2]) %>%
layout(title="Valores ajustadosvs Valores observados")
p
No trace type specified:
Based on info supplied, a 'scatter' trace seems appropriate.
Read more about this trace type -> https://plot.ly/r/reference/#scatter
No scatter mode specifed:
Setting the mode to markers
Read more about this attribute -> https://plot.ly/r/reference/#scatter-mode
No trace type specified:
Based on info supplied, a 'scatter' trace seems appropriate.
Read more about this trace type -> https://plot.ly/r/reference/#scatter
No scatter mode specifed:
Setting the mode to markers
Read more about this attribute -> https://plot.ly/r/reference/#scatter-mode
Podemos ver que como el modelo tiene una \(R^2 ~ .8592\) por lo que podemos decir que el modelo esta bien especificao ya que explica el aproximadamente el 85% de la varianza con las variables que estamos usando. nuestro error estandar residual \(\sigma=1497\) por lo que \(sigma^2=1497^2\)
\[arcos(R^2)=\theta\] En este caso en particular \(R^2=.8592\) por lo que
acos(.8592)
[1] 0.5370923
En caso de que podamos aproximar con una normal tenemos lo siguiente:
mu<-X%*%beta
Error in X %*% beta : non-conformable arguments
En nuestro caso en particular:
X_reg<-cbind(rep(1,length(y)), diamonds_num[,2:7])
y_reg<-diamonds_num[,1]
betasigma<-c(0,0,0,0,0,0,0,80)
optim(par=betasigma,lik1, X=X_reg, y=y_reg)
Error in optim(par = betasigma, lik1, X = X_reg, y = y_reg) :
objective function in optim evaluates to length 0 not 1